/*
 * Preprocessor2D.cpp
 *
 *  Created on: Aug 26, 2011
 *      Author: tpz22rc
 */

#include "Preprocessor2D.h"
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <cmath>

using std::ifstream;
using std::clog;
using std::cerr;

Preprocessor2D::Preprocessor2D() {
	numOfC=numOfC=numOfB=0;
}

Preprocessor2D::~Preprocessor2D() {
}

Preprocessor2D::Preprocessor2D(const int &numOfC,const int &numOfE,const int &numOfB){
	this->numOfC=numOfC;
	this->numOfE=numOfE;
	this->numOfB=numOfB;
}

double ** Preprocessor2D::readCoord(const char *f){
	double **coord=new double*[numOfC+1];
	for(int i=0;i<=numOfC;i++)
		coord[i]=new double[MAX_COORD+1];
	ifstream in(f);
	if(!in.is_open()){
		std::cerr<<"error reading file: "<<f<<"\n";
		abort();
	}
	for(int i=1;i<=numOfC;i++){
		for(int j=1;j<=MAX_COORD;j++){
			in>>coord[i][j];
		}
	}
	in.close();
	return coord;
}

int ** Preprocessor2D::readElement(const char *f){
	int **connect=new int*[numOfE+1];
	for(int i=0;i<=numOfE;i++){
		connect[i]=new int[MAX_CONNECT+1];
	}
	ifstream in(f);
	if(!in.is_open()){
		std::cerr<<"error reading file: "<<f<<"\n";
		abort();
	}
	for(int i=1;i<=numOfE;i++){
		for(int j=1;j<=MAX_CONNECT;j++){
			in>>connect[i][j];
		}
	}
	in.close();
	return connect;
}

int ** Preprocessor2D::readBound(const char *f){
	int **bound=new int*[numOfB+1];
	for(int i=0;i<=numOfB;i++)
		bound[i]=new int[MAX_BOUND+1];
	ifstream in(f);
	if(!in.is_open()){
		std::cerr<<"error reading file: "<<f<<"\n";
		abort();
	}
	for(int i=1;i<=numOfB;i++){
		for(int j=1;j<=MAX_BOUND;j++){
			in>>bound[i][j];
		}
	}
	in.close();
	return bound;

}

double ** Preprocessor2D::geometry(double **coord,int **connect){
	double x[MAX_CONNECT + 1], y[MAX_CONNECT + 1];
	double area;
	int cp;

	double **geo = new double*[numOfE + 1];
	for (int i = 0; i <= numOfE; i++)
		geo[i] = new double[GEO + 1];

	for (int e = 1; e <= numOfE; e++) {
		for (int j = 1; j <= MAX_CONNECT; j++) {
			cp = connect[e][j];
			x[j] = coord[cp][1];
			y[j] = coord[cp][2];
		}
		geo[e][1] = y[2] - y[3];
		geo[e][2] = y[3] - y[1];
		geo[e][3] = y[1] - y[2];
		geo[e][4] = x[3] - x[2];
		geo[e][5] = x[1] - x[3];
		geo[e][6] = x[2] - x[1];
		area =(x[2]*y[3] - x[3]*y[2]) + (x[3]*y[1] - x[1]*y[3])+ (x[1]*y[2] - x[2]*y[1]);
		geo[e][7] = area / 2;

	}
	return geo;
}

double ** Preprocessor2D::getnorm(double **coord,int **bound){
	double **norm = new double*[numOfB + 1];
	int j;
	for (j = 0; j <= numOfB; j++)
		norm[j] = new double[MAX_BOUND + 1];
	int i, ip0, ip1;
	double dx, dy, rl;
	for (i = 1; i <= numOfB; i++) {
		ip0 = bound[i][1];
		ip1 = bound[i][2];
		dx = coord[ip1][1] - coord[ip0][1];
		dy = coord[ip1][2] - coord[ip0][2];
		rl = sqrt(dx * dx + dy * dy);
		norm[i][1] = dy / rl;
		norm[i][2] = -dx / rl;
		norm[i][3] = rl;

	}
	return norm;
}

double * Preprocessor2D::lmass(double **geo,int **connect){
	double *lm = new double[numOfC + 1];
	int i;
	for (i = 1; i <= numOfC; i++)
		lm[i] = 0.0;
	int ielem, inode, in;
	double rj, rj6;
	for (ielem = 1; ielem <= numOfE; ielem++) {
		rj = geo[ielem][GEO];
		rj6 = rj / 3;
		for (inode = 1; inode <= MAX_CONNECT; inode++) {
			in = connect[ielem][inode];
			lm[in] = lm[in] + rj6;
		}
		for (inode = 1; inode <= numOfC; inode++) {
			lm[inode] = 1.00 / lm[inode];
		}

	}
	return lm;
}
